%% Organize DATA

names = fullnames_data(choose_variab)
bettername = betternames_data(choose_variab) 
region = region_data(choose_variab)
names_exog = fullnames_data(choose_exog)
Y = data_all(:,choose_variab);
Y_exog=data_all(:,choose_exog);


% chose names for IRFs
names_IRFs = fullnames_data(choose_variab);
transformIRF_used = transformIRFs_data(choose_variab);
labelIRF_used = labelIRFs_data(choose_variab);
Numbers = 1:length(day);
Numbers=Numbers';

% reduce to period for which all measures are available
cond = isnan(Y);
cond2 = sum(cond,2) == 0;
Y = Y(cond2,:);
Y_exog = Y_exog(cond2,:);
year_available     = year(cond2);
month_available    = month(cond2);
day_available    = day(cond2);
Numbers_available = Numbers(cond2);

data_all_available = data_all(cond2,:);
% keep only after a certain period, in case data are available for all variables added for before that period

% full sample
from = 1 
until = 30000; 

cond = and(Numbers_available >= from,Numbers_available <= until);
sum(cond)
Y = Y(cond,:);
Y_exog = Y_exog(cond,:);

year_available = year_available(cond);
month_available = month_available(cond);
day_available    = day_available(cond);
Numbers_available = Numbers_available(cond);

disp(['-->> Data available for :    ', num2str(year_available(1)) 'M' num2str(month_available(1)) 'D' num2str(day_available(1)) ...
                            ' until ', num2str(year_available(end)) 'M' num2str(month_available(end)) 'D' num2str(day_available(end))] );
% ensure Y enters horizontally
Y = Y';
Y_exog = Y_exog';
k = size(Y,1)
T = size(Y,2)
assert(  T > k)
assert( length(year_available) == T)      

year_available_help=zeros(T,1);
for i=1:length(year_available)
year_available_help(i,:) = year_available(i,:) + 0.001*i;
end

%% Estimate RVAR
maxlags = 20;
 
constantD   = 1;
trendD      = 0;
trendBreaks = [];
seasonD     = 0;
moreBreaks  = [];

[Ahat, Resid, Sigmahat, Ahat_determ, X_determ, T_used, df_lost] = ...
         zFC_VAR_Estimation_exog(Y,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);        
        
year_available_used     = year_available(p+1:end);
assert(length(year_available_used) == size(Resid,2))
month_available_used    = month_available(p+1:end);
day_available_used    = day_available(p+1:end);
Numbers_available_used = Numbers_available(p+1:end);

%% IRFs, point estimation
T_step = [1:1:T_irf]; % used to plot IRFs
shock_considered = 1; % variable shocked

% CHOLESKY
if cholesky_on == 1
    B_chol          = chol(Sigmahat)';
    shock_size_chol = 1/B_chol(1,1);
    IRF_chol        = zFC_IRFs(Ahat, B_chol, T_irf, shock_considered, shock_size_chol);
end 

Instrument = data_all(:,choose_shock);
cond_i = isnan(Instrument);
cond2_i = sum(cond_i,2) == 0;
Instrument = Instrument(cond2_i,:);

year_available_ins     = year(cond2_i);
month_available_ins    = month(cond2_i);
day_available_ins = day(cond2_i);
Numbers_available_ins = Numbers(cond2_i);

% select subperiod common to shocks used and variables used
from2_year         = max(year_available_used(1),year_available_ins(1));
from2_month          = max(month_available_used(1),month_available_ins(1));
from2_day          = max(day_available_used(1),day_available_ins(1));
from2          = max(Numbers_available_used(1),Numbers_available_ins(1));
% % % % %from2          = timeaxis_used(59);
until2_year         = min(year_available_used(end),year_available_ins(end));
until2_month         = min(month_available_used(end),month_available_ins(end));
until2_day         = min(day_available_used(end),day_available_ins(end));
until2         = min(Numbers_available_used(end),Numbers_available_ins(end));
Numers_available_used2 = [from2:until2];
%                         
disp(['-->> Period effectively used in the identification:    ', num2str(year_available_ins(1)) 'M' num2str(month_available_ins(1)) 'D' num2str(day_available_ins(1)) ...
                            ' until ', num2str(year_available_ins(end)) 'M' num2str(month_available_ins(end)) 'D' num2str(day_available_ins(end))] );
                      
% Take only the instruments in this subperiod
cond1            = Numbers_available_ins >= from2 & Numbers_available_ins <= until2;
Instruments_used = Instrument(cond1,:);
n_instrum=size(Instrument);
n_instrum=n_instrum(1,2);
 

% Take only the residuals in this subperiod
cond2      = Numbers_available_used >= from2 & Numbers_available_used <= until2;
Resid_used = Resid(:,cond2);
AAA=[Resid_used' Instruments_used];

names_instrument = fullnames_data(choose_shock) ;
[bla,g] = size(Instruments_used) ;

% Identify the model
% choose correct position for indicator variable
position_indicator = 1;

[imp_vect_relat, absolute_impulse_vector, ...
             F_stat_test1, p_values_test1, signif_test1, Rsquared_test1, Betas_test1, ...
             F_stat_test2, p_values_test2, signif_test2, Rsquared_test2, ...
             correl_instrument_shock, a_vect, estim_shocks] = zFC_IdentifExtInstruments5(Resid_used, Instruments_used, position_indicator, Sigmahat);
% 
table_validity_rows = char('beta', 'numb. of stars', 'numb. of obs.', 'F statistic', 'R squared', 'obs. of the proxy used');
table_validity_rows = cellstr(table_validity_rows); 
table_validity = [Betas_test1'; signif_test1(:,1)'; length(Instruments_used)*ones(1,k)/1000; F_stat_test1'; Rsquared_test1']    

nonzero_instrument = nonzeros(Instruments_used);
n_nonzero_instrument = length(nonzero_instrument)

B_extin    = [absolute_impulse_vector, zeros(k,k-1)];

if sd_shock == 1
    shock_size_extin = 1;
else
    shock_size_extin = shock_size/absolute_impulse_vector(1);
end
        
IRF_extin  = zFC_IRFs(Ahat, B_extin, T_irf, 1, shock_size_extin);


%% Bootstrap on the reduced form model 

if bootstrap_on == 1 
    
    Yinit           = Y(:,1:p);
    Ahat_loop       = NaN*ones(k,k*p,Loops);
    Resid_used_loop = NaN*ones(k,length(Resid_used),Loops);
    Sigmahat_loop   = NaN*ones(k,k,Loops);
    Instruments_used_loop = NaN*ones(length(Instruments_used),n_instrum,Loops);

    for l = 1:Loops

        % Generate a vector that flips signs randomly
        step = randn(T_used,1);
        change_sign = -1*ones(T_used,1);
        change_sign(step > 0) = 1;

        % Generate pseudo data
        Resid_pseudo = Resid*diag(change_sign);

        Y_pseudo = zFC_GenerateVARData(Yinit, Ahat, Resid_pseudo, Ahat_determ, X_determ);

        % Estimate the VAR on the pseudo data
        [Ahat_loop(:,:,l), Resid_step, Sigmahat_loop(:,:,l), ~, ~, ~, ~] = ...
                       zFC_VAR_Estimation_exog(Y_pseudo,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);

        % Store residuals for the subpart that overlaps with the shocks used
        Resid_used_loop(:,:,l) = Resid_step(:,cond2);

        % Flip sign of the instruments using the same order in the bootstrapped residuals, consider only the part that overlaps with the residuals
        Instruments_used_loop(:,:,l) = diag(change_sign(cond2))*Instruments_used;
    end

    % IRFs, error bands
    IRF_chol_loop  = NaN*ones(k,T_irf,Loops);
    IRF_extin_loop = NaN*ones(k,T_irf,Loops);
    IRF_extin_new_loop = NaN*ones(k,T_irf,Loops);

    for l = 1:Loops

        if cholesky_on == 1
        % Cholesky
        B_chol_loop = chol(Sigmahat_loop(:,:,l))';
        shock_size_chol_loop = shock_size_chol;
        IRF_chol_loop(:,:,l)  = zFC_IRFs(Ahat_loop(:,:,l), B_chol_loop, T_irf, shock_considered, shock_size_chol_loop);
        end

        % External instrument, using only monthly data
        [~, ~, relative_impulse_vector_loop, ~, incons_impulse_vector_loop] = zFC_IdentifExtInstruments2(Resid_used_loop(:,:,l),Instruments_used_loop(:,:,l));
        b11_loop = zFC_ScaleInExternalInstrument(Sigmahat_loop(:,:,l), relative_impulse_vector_loop);
        absolute_impulse_vector_loop = relative_impulse_vector_loop*b11_loop;

        B_extin_loop    = [absolute_impulse_vector_loop, zeros(k,k-1)];
        
        if sd_shock == 1
            shock_size_extin_loop = 1; 
        else
         shock_size_extin_loop = shock_size/absolute_impulse_vector_loop(1);
        end

        IRF_extin_loop(:,:,l)  = zFC_IRFs(Ahat_loop(:,:,l), B_extin_loop, T_irf, 1, shock_size_extin_loop);

    end

    % Compute error bands
    if cholesky_on ==1
        IRF_chol_bands = NaN*ones(k,T_irf,2);
    end

    IRF_extin_bands = NaN*ones(k,T_irf,2);
    interval       = 10;
    for ttt = 1:T_irf 
        for kkk = 1:k
            if cholesky_on == 1
                IRF_chol_bands(kkk,ttt,1)   =  prctile(IRF_chol_loop(kkk,ttt,:),interval/2);
                IRF_chol_bands(kkk,ttt,2)   =  prctile(IRF_chol_loop(kkk,ttt,:),100-interval/2);   
            end
            IRF_extin_bands(kkk,ttt,1)  =  prctile(IRF_extin_loop(kkk,ttt,:),interval/2);
            IRF_extin_bands(kkk,ttt,2)  =  prctile(IRF_extin_loop(kkk,ttt,:),100-interval/2);   
        end
    end

end

%% model with cumulated first day of intervention and Cholesky deomposition

choose_variab = [13,2:8]; % Cholesky with cumulated first days of sequencs 

names = fullnames_data(choose_variab)
bettername = betternames_data(choose_variab) 
region = region_data(choose_variab)
names_exog = fullnames_data(choose_exog)
Y = data_all(:,choose_variab);
Y_exog=data_all(:,choose_exog);


% chose names for IRFs
names_IRFs = fullnames_data(choose_variab);
transformIRF_used = transformIRFs_data(choose_variab);
labelIRF_used = labelIRFs_data(choose_variab);
Numbers = 1:length(day);
Numbers=Numbers';

% reduce to period for which all measures are available
cond = isnan(Y);
cond2 = sum(cond,2) == 0;
Y = Y(cond2,:);
Y_exog = Y_exog(cond2,:);
year_available     = year(cond2);
month_available    = month(cond2);
day_available    = day(cond2);
Numbers_available = Numbers(cond2);
% plot(Y(:,2));

% % % comparison with Cholesky using only first day of intervention to
% % % construct cumulated reserves
Y(:,1) = cumsum(Y(:,1)) ;
% % % 

data_all_available = data_all(cond2,:);
% keep only after a certain period, in case data are available for all variables added for before that period

% full sample
from = 1 % + 2610 + 00;
until = 30000; % end 2015

cond = and(Numbers_available >= from,Numbers_available <= until);
sum(cond)
Y = Y(cond,:);
Y_exog = Y_exog(cond,:);

year_available = year_available(cond);
month_available = month_available(cond);
day_available    = day_available(cond);
Numbers_available = Numbers_available(cond);

disp(['-->> Data available for :    ', num2str(year_available(1)) 'M' num2str(month_available(1)) 'D' num2str(day_available(1)) ...
                            ' until ', num2str(year_available(end)) 'M' num2str(month_available(end)) 'D' num2str(day_available(end))] );
% ensure Y enters horizontally
Y = Y';
Y_exog = Y_exog';
k = size(Y,1)
T = size(Y,2)
assert(  T > k)
assert( length(year_available) == T)      

year_available_help=zeros(T,1);
for i=1:length(year_available)
year_available_help(i,:) = year_available(i,:) + 0.001*i;
end


% Estimate RVAR
maxlags = 20;
 
constantD   = 1;
trendD      = 0;
trendBreaks = [];
seasonD     = 0;
moreBreaks  = [];

[Ahat, Resid, Sigmahat, Ahat_determ, X_determ, T_used, df_lost] = ...
         zFC_VAR_Estimation_exog(Y,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);        
        
year_available_used     = year_available(p+1:end);
assert(length(year_available_used) == size(Resid,2))
month_available_used    = month_available(p+1:end);
day_available_used    = day_available(p+1:end);
Numbers_available_used = Numbers_available(p+1:end);

% IRFs, point estimation
T_step = [1:1:T_irf]; % used to plot IRFs
shock_considered = 1; % variable shocked

% CHOLESKY
if cholesky_on == 1
    B_chol          = chol(Sigmahat)';
    shock_size_chol = 1/B_chol(1,1) ;
    IRF_chol2        = zFC_IRFs(Ahat, B_chol, T_irf, shock_considered, shock_size_chol);
end 

% Identify the model
% choose correct position for indicator variable
position_indicator = 1;

%% Final figure

figure(7) 

for i = 1:k
    subplot(rows_fig,cols_fig,i)

    if bootstrap_on == 1
    
    fill([T_step,fliplr(T_step)],[transformIRF_used(i)*IRF_extin_bands(i,:,2), fliplr(transformIRF_used(i)*IRF_extin_bands(i,:,1))], ...
                        grey,'EdgeColor','none', 'FaceAlpha', 0.5), hold on
    end
    
     a = plot(T_step, transformIRF_used(i)*IRF_extin(i,:),     'k-', 'Linewidth', 2);  hold on

     if cholesky_on == 1
         fill([T_step,fliplr(T_step)],[transformIRF_used(i)*IRF_chol_bands(i,:,2), fliplr(transformIRF_used(i)*IRF_chol_bands(i,:,1))], ...
                            blue,'EdgeColor','none', 'FaceAlpha', 0.5), hold on
          b = plot(T_step, transformIRF_used(i)*IRF_chol(i,:),     'b--', 'Linewidth', 2) ; hold on
          

     end
    title_i = strcat(bettername(i)) ;
    plot(T_step,zeros(T_irf,1), 'k')
    title(title_i,'FontSize',font_size)
       ylabel(labelIRF_used(i),'FontSize',font_size)
    axis tight
    set(gca,'box','off','FontSize',font_size)
    set(gca,'box','off')
    if cholesky_on == 1
        if i == 1
           legend([a b], 'Signaling and portfolio','Portfolio channel', 'Location', 'Southeast'), legend('boxoff') 
        end
    end
end

